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ABSTRACT 

We present novel statistical tools to cross-correlate frequency cleaned thermal Sunyaev-Zel'dovich 
(tSZ) maps and tomographic weak lensing (wl) convergence maps. Moving beyond the lowest or- 
der cross-correlation, we introduce a hierarchy of mixed higher-order statistics, the cumulants and 
cumulant correlators, to analyze non-Gaussianity in real space, as well as corresponding polyspec- 
tra in the harmonic domain. Using these moments, we derive analytical expressions for the joint 
two-point probability distribution function (2PDF) for smoothed tSZ (ys) and convergence (k^) 
maps. The presence of tomographic information allows us to study the evolution of higher order 
mixed tSZ-weak lensing statistics with redshift. We express the joint PDFs pKy{Ks, Us) in terms of 
individual one-point PDFs (pk(ks), Py{ys)) and the relevant bias functions (6„(ks), by{ys)). Ana- 
lytical results for two different regimes are presented that correspond to the small and large angular 
smoothing scales. Results are also obtained for corresponding hot spots in the tSZ and convergence 
maps. In addition to results based on hierarchical techniques and perturbative methods, we present 
results of calculations based on the lognormal approximation. The analytical expressions derived 
here are generic and applicable to cross-correlation studies of arbitrary tracers of large scale struc- 
ture including e.g. that of tSZ and soft X-ray background. 
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1 INTRODUCTION 

Free electrons in the Universe can be detected through the inverse Compton scattering of Cosmic Microwave Background (CMB) photons 
dSunvaev & Zeldovich|[l980ll 19721; lBirkinshawlll999l ; lReDhaelilll995h . An inhomogeneous distribution of electrons thus induces a secondary 
anisotropy in the CMB radiation, which is proportional to the thermal energy of the electrons integrated along the line of sight directio n. This well- 
known effect is called the thermal Sunya ev-Zel'dovich effect (tSZ) and it is now routinely being used to image nearby galaxy clusters ( iReese et"al] 
I2OO2I : Ijones et al.ll2005l : iLaRogudllOOa) which the electron temperature can reach the order of lOkeV. The future seems set for a huge increase 
in Sunyaev-Zel 'dovich measurements of ga laxy clusters; early data from Planck alone has generated 189 candidates with reasonably high signal- 
to-noise ratios jPlanck Collaborationll201ll) . In the intergalactic medium (IGM), however, the gas is expected to be in mildly overdense regions 
and to reach a temperature of only IkeV or so. While the tSZ effect from clusters can alter the temperature of the CMB by an amount of order 
of mK in the Rayleigh- Jeans part of the spectrum, the contribution from the IGM is expected only to reach the /iK range. This is below the 
threshold for detection by WMAP, but sho uld nevertheless be detectable by the on going Planck Survey oi0 due to its wider frequency coverage, 
higher sensitivity, and improved resolution l lHansen et all2005l : lloudaki et all2010l) . 

The unique spectral dependence of the tSZ effect helps in the task of separating it from other sources of CM B temperature fluctuations so var- 
ious well developed component separation schemes are available for construction of frequency-cleaned SZ maps l lLeaclj2008l ; lBouchet & Gisperj 
1 1999 : Delabroullie, Card oso & Patanchon 2003). The construction of such maps will provide us with a direct opportunity to probe the thermal 
history of the Universe, in tandem with other obse rvations: owing to their thermal ( peculiar motions) ionized electrons scatter CMB photons, 
an effect which can be studied using th e tSZ effect (iRosati. Borgani & Normanll 19981) : the neutral component of the IGM can be also be probed 
via observations of the Lyman-a forest l lRauch|[T998l) : X-ray emission due to the thermal bremmstrahlung that originates from the interaction of 
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ionized electrons and nuclei can also provide clues to the nature and evolution of the IGM. These tracers all probe different phases of the IGM and 
thus play complementary roles. For example, the X-ray emission depends on the square of the density, so is most sensitive to the densest regions 
in the IGM, primarily in the local Universe, while tSZ studies can probe the more distant Universe because Compton scattering is independent of 
redshift, as well as being an unbiased tracer of all electrons, as they all participate in the scattering, the tSZ effect is weakly dependent on density 

so probes the electron density in a wide range of environments. 

The modeling of lower order statistics of the tSZ effect can be carri ed out using various approaches. In the past Seliakll2000l: Zhang & PenI 
l200ll ; lKomatsu & SeliaklboOlUZhang & Setri[2007l : ICooravll2000l . l200lh the modeling has followed the halo model Jcoorav & Sethll2002ll . The 
statistical distribution of haloes, specificall y their number-density as a function of mass (i.e. their mass function), is assumed to be that predicted 
provid ed by the Press-Schechter forma lism jPress & Schechteill974h or its generalisations, and the radial profile of such haloes were assumed to be 
that o f iNavarro. Frenk & White! h996l) . The hot gas, assumed to have been heated by shocks, is taken to be in hydrodynamical equilibrium with the 
dark matter distribution. The typical temperature reached in such syste ms is sufficien t to ionize the hydrogen and helium atoms. These ingredients 
are sufficient to model the tSZ effect raising from collapsed haloes l lCooravlboOOl) ; in addition to thi s analytical modeling, numerical Simula 



tion of the tS Z plays an important role in our understa nding of th e physics involved dPersi et al 
Springeletafc oOl: White, Hernquist & Springel 2002: Lin et alj |2004: Zha ng, Pen & Tracll2004 



ll995l:lRefre gier et al."2000':'Seliak et al."200I 



Cao. Liu & Fang 2007; Roncarelli et al. 2007 



Hallman et al.ll2009l , 120071 : Ida Silva et al.ll2000l) . To extend the halo model to larger scales, the extended distribution of free electrons is typically 



assumed to trace the distribution of dark matter on large scales where the variations in density are in the linear or quasi-linear regime. A perturba- 
tive ap pr oach along these line s has been developed by several authors over the years ( ICooravl2000ll200ll : ICoorav & Hul2000l : lGoldberg & Spergell 
ll999allbllMunshi et al.l20lTal lbl). 

Ongoing and proposed ground surveys su ch as sz^E actE ape}<3 sptQ and the ongoing all sky survey Planck will map the y-sky with 
a gre at precision. The SPT ([Lueker et al.ll2010h collaboration has already reported the measurement of the tSZ power spectrum at I ~ 3000; the 
ACT jFowler et alJl2010l ; lDunklev et al.ll2010^ collaboration has also reported analysis on similar scales. It is expected that ongoing surveys will 
improve these measurements due to their improved sky coverage as well as wider frequency range. It is important to realize why the study of 
secondaries such as tSZ should be an important aspect of any CMB mission. In addition to the important physics the secondaries probe, accurate 
modeling of the secondary no n-Gaussianities is required to avoid 20% — 30% constraint degradations in future CMB data-sets such as Planck and 
CMBPofl dSmidt et allioToh . 

While the tSZ surveys described above provide a direct probe of the baryonic Universe, weak lensing observations jMunshi et al.ll200^ on 
the other hand can map the dark matter distribution in an unbiased way. In recent years there has been tremendous progress on the technical 
front in terms of specification and control of systematics in weak lensing observables. There are many ongoing and future weak lensing surveys 

such as cfhtQ legacy survey. Pan- STARR^ and the Dark Energy survey (DESfl and in future, the Large Synoptic Survey Telescope (LSSTj |, 

Joint Dark Energy Mission (jdem]3 will map the dark matter and dark energy distribution of the entire sky in unprece dented details. Among 
other things, hold great promise in shedding light on the nature of the dark energy and the origin of neut rino masses I Joudaki & Kaplinghall 
I20T J), where the weak lensing signals dominate the others considered by e.g. the Dark Energy Task Force jAlbrech et alj|20l"]|). However, the 
optimism associated with lensing is predicated on overcoming the va st systematic uncertai nties in both measurement and in theo ry ( Hirata & Se lia^ 
I2OO4I ; iMa. Hu & Huteredl2005l ; Tcoorav & Hull2002l ; iHirata & Seli ak 2003: White., 20041 : iHuteredlloO^ : [McDonald. Trac & Contaldi^ 2006). The 
statistics of the weak lensing conv ergence have been studied in great det ail using an extension of perturbation theory |Munshi & Jainll2000ll200ll : 
iMunshi. Valageas & Barbeil2004l) and methods based on the halo model jCoorav & Hij200"ol : lTakada & Jainl2002l . l2003 ). These studies developed 
techniques that can be used to predict the lower-order moments (equivalent to the power spectrum and multi-spectra in the harmonic domain) and 
the entire PDF for a given weak lensing survey. The photometric redshifts of source galaxies a re useful for tomographic studies of the dark matter 
distribution and establish a three-dimensional picture of their distribution jMunshi et alj2OI0l) . 

Due to the line of sight integration inherent in the tSZ effect, the redshift evolution is completely lost which degrades the power of tSZ to 
distinguish different thermal histories. It is however possible to recover some of the information lost by cross-correlating the tSZ with external 
tracers (Shao et al. 201 la; Zhang & Pen 2001). These tracers could comprise galaxy surveys with spectroscopic redshift information or dark matter 
distribution from weak lensing surveys in tomographic slices. Such cross-correlations has been studied using two-point sta tistics. We extend these 
results to higher order. Higher order statistics such as cumulants and cumulant correlators jMunshi. Melott & Coles|[r999l) can probe higher order 
cross-correlations and can in principle provide an independent probe of the bias associated with the baryonic pressure fluctuations. The tomographic 
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Figure 1. The parameter Km, that denotes the minimum value of convergence is displayed as a function of source redshift for the ACDM cosmology (right panel). 
The parameter Km is independent of angular smoothing scale. The dashed line in the right panel corresponds to Km for the projected 2D survey. The left panel 
shows the source distribution as a function of redshift. 



cross-correlation statistics that we develop here can be applied to surveys with overlapping sky coverage. Many of the weak lensing surveys and the 
tSZ surveys will have overlapping sky coverage. For example, DBS will have overlap with the SPT sky coverage and plans to measure photometric 
redshifts of roughly 10** galaxies up to z = 1.3. Tomographic cross-correlation statistics at second or higher-order can provide a more detailed 
picture of the evolution of the thermal history of the baryonic gas by mapping the associated pressure fluctuation. The results presented here are 
generic and can be extended to study other secondaries, such as the cross-correlation involving CMB lensing and tSZ maps. 

This paper is organized as follows. In Sj2]we introduce our notations for both tSZ and weak-lensing observables. Next we introduce the mixed 
cumulant correlato rs and the cumulants in i|3] Two different models are used to model the underlying dark matter distribution. We employ the 
dSmith et al.ll2003l) prescription, for modelling the evolution of matter power spectrum and use the hierarchical approach for modelling of the 
higher order correlation hierarchy. This allows us to extend the lower order moment results to arbitrary order and construct the relevant PDF and 
bias for the smoothed weak lensing k field and the y maps in S|4] A brief review of the hierarchical ansatz in the quasilinear and non-linear regimes, 
as well as the lognormal approach, is provided in the appendix. 



2 NOTATIONS 

In this section, we introduce our notation for the tSZ effect and weak lensing convergence. We will use the following line element: 

ds^ = -c^dt^ + a^(r) [dr^ + d\{r){sm OdO"^ + d<p'^)] (1) 

Here dA{r) is the comoving angular diameter distance at a (comoving) radial distance r. The angular diameter distance is dA{r) = 
(— K)"^''^ sinh((— K)^''^r) for negative curvature, dA{r) = (K)"^'''^ sin((K)^''^r) for positive curvature, and for a flat Universe dA{r) = r. 
The parameters Ho and Oo decide the constant K: K = (Slo — l)Ho. The underlying cosmology that we adopt for numerical study is specified 
by the following parameter values (to be introduced later): SIa ~ 0.741, h = 0.72, ilb = 0.044, ficDM = 0.215, f2o = f^ii + f^CDM, ?^s = 
0.964, Wo ~ —1, Wa = 0, as = 0.803, = 0. The comoving radial distance at a redshift z is determined by the cosmology (ilo. Ho) 

r(z)= "'^ (2) 

Jo Hov/f^o(l + + nK{l + zY + S^A 

Throughout, c will denote the speed of light and will be set to unity. 



2.1 Thermal Sunyaev Zel'dovich Effect 

The tSZ effect contributes to the CMB temperature fluctuation and is typically expressed as iriy,^) = Ar(f2)/To = g{x^)y{Q,). In this 
expression, g{xi,) corresponds to the spectral dependence and encodes the angular dependence; Xi, represents the dimensionless frequency 
and Q. — (6, (f>) corresponds to a unit vector that signifies pixel positions. A subscript s will be used to denote the smoothed maps, e.g. ysiCl). In 
the non-relativistic limit g{xv) takes the following form: 



g{x,,) = x^ coth (^y^ - 4 = {^^ ^^^^ ^ I - 4^ 



hv V 5.28mm 

fcsTo " 56.84GHz " A ' ^ ' 
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Here and h are the Boltzmann and Planck constants respectively; u denotes the frequency of the photon and To = 2.726 K is the temperature 
of the CMB sky. The tSZ effect shows as CMB temperature decrement at v 218 GHz and as an increment at i/ ^ 218 GHz, with a null at 
v = 218 GHz. In the Rayleigh Jeans limit, characterized by x <^ 1, g(x) ~ —2 is roughly independent of frequency and in the other limiting 
situation Xi, ^ 1, g{x,^) ~ (x^ — 4). The key information regarding thermal history of the Universe is encoded in y{Cl) maps that are extracted 
from frequency maps obtained through multifrequency CMB observations. The y maps are opacity weighted integrated pressure fluctuations along 
the line of sight. 

yiyi) = ds UeCTT 7 = 7 / dr a{r)nekBTe{^,r) = / dr) ahj) TlJ-q^Q.) = / dr wszMneir). (4) 

J nieC^ nieC^ Ja m^c^ Jo Jo 

Here He = UeksTe, while rrie conesponds to the electron mass, ks denotes the Boltzmann's constant, ar ~ 6.6510^^^cm^ represents the 
Thompson cross-section, denotes the number-density of electrons, and Te is the electron temperature. Conformal time is denoted by dr] = 
dt/a{t). The line-of-sight integral depends on the comoving radial coordinate distance r and a(r-) is the corresponding scale factor of the Universe. 
The weight is defined as wsz{r) = T(r) — aTne{r)a{r), where the dot defines the derivative with respect to comoving radial distance r, and the 
3D pressure fluctuation is defined as tte = kBTe/meC^. We will cross-correlate the comptonization map with tomographic and projected 
maps from weak lensing surveys to constrain the thermal history of the Universe and its evolution with redshift. Throughout, we will consider 
the Rayleigh-Jeans part of the spectrum St = — 2y; for ACT and SPT operating at v — 150GHz from Eq. ^ we get g{x) — —0.95. Detailed 
modelling of the bias is required only for the computation of variance. The variance {Sy'^{Cl)) samples the pressure fluctuation power spectrum 
Ptttt and is expressed as: 

The pressure power spectrum P,rir {k, z) at a redshift z is expressed in terms of the underlying power spectrum P^^ (k, z) using a bias fe^ (k, z), i.e. 
P,r7r(fc, z) = z)¥si{k, z). The bias i),r(fc, z) is assumed to be independent of length scale or equivalently wave number k, i.e. 6^(A;, z) = 
hTT{z). The redshift dependent bias can be expressed as b,r (z) = bir (0)/(l + z). Here 6,r (0) can be written as fe,r (0) = kBTe(Q)bs / m e (? ■ D if feren t 
values of 6^ have been reported, e.g. feefregier et al. 2000 ) found &j ^ 8 - 9 and Te(0) ^ 0.3 - 0.4. On the other hand (Selj ak et al, ; 200lh 
found hs-i-A. Typical v alue of fO) found by llCen & OstrikeJI 1 9991) is 6^(0) = 0.0039. This is a factor of two lower than the value used by 
l lGoldberg & SDergelll l999al lbl) and l lCoorav & Htj|2000l) . A Gaussian beam bi{6s) with FWHM at 9s is assumed. The PDF and bias functions of 
hot-spots are specified in terms of the normalized cumulants or cumulant correlators, which we probe next. 



(5) 



2.2 Weak Lensing in Projection and Tomography 

Cosmological weak lensing effects are conveniently encoded in the effective convergence field, which is defined as a weighted projection of the 
matter density contrast S. The statistics of the smoothed weak lensing convergence Ks are similar to the projected 3D density contrast (5(x) and 
can be expressed through a line of sight integration — J^" dr 'w^i{r)S{r, Cl). Their weight w^i{r) is sensitive to the source distribution; 

for 2D surveys without any source redshift information, can be expressed as 

w^i{r) = ^^Q.Ma~^{T)dA{r) j dr ps(r) ^ '^^^''^^ \ Psizs) ^ exp'''^""" . (6) 
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Figure 3. The redshift-resolved cross spectra (7''''°'^ given in Eq. I12t for weak lensing tomographic slices and the tSZ survey are plotted as function of harmonics 



Here P3{zs) represents the source redshift distribution. The peak of the distribution is reached at 2zo. We will adopt two different survey config- 
urations. We adopt zq — 0.5 for the LSST survey. For the tomographic calculation we divide the entire source redshift range in five redshift bins 
with each tomographic bin containing roughly the same number density of source galaxies. The redshift bins of the sources are delimited at [0.75, 
1.1, 1.45, 1.95, 3.00]. The tomographic convergence maps k'*' — JJ^° dr w^' {r)S{r, Cl) depend on individual weights: 



^ilM = ?^f^M-^a-^(r) dA{r) 



2 c2 



Hi 



dr ps{r] 



dzs dA{r — rs) 



dr. 



dij 



(7) 



The projected convergences for individual tomographic bins are defined as = — j^' dr w^^l (r). We will next use these expressions to derive 

the cross-correlations at second- and higher-order both in real space and the harmonic domain. 

Some comments a r e in o r der at th i s poin t . The hierarchical model that we use here was previously used by IZhang. Pen & Wan3 1 I2OO2I) : 
ICoorav. Tegmark"& Hul |2003); ICooravl l l200d. 1200 ih . It is known that the j/s -parameter is proportional to the square of the density contrast 
5^, so e.g. the three- point correlation func tion of i/s effectively samples the six-point correlation function of the underlying density contrast 5. 
However, studies by l lCooravll2000l,l200lh assumed a linear biasing model for the pressure fluctuations, i.e. in the Fouri er domain 5Tr(k,z) = 
5{k, z)bTr{k, z). The pressure bias model assumes the pressure to be a linear convolution over the density field. Finally ICooravl l l200(3l l200lh 
further simplified the bias 6-n-(fc, z) by factorizing it into a redshift dependent and momentum k dependent part. The spatial dependence was 
assumed to be constant. This is the simplification that we will use in our results, which will be useful in computation of correlatio n hierarchy to 
an arbitrary order. The hierarchica l ansatz and the resulting scaling functions were also used by jValageas. Schaefer & Silk 2001,; Valageas & Silkl 
ll999l : IValageas & Schaeffej2000h . For computation of the number density of collapsed objects these were then used to obtain the thermal Sunyaev- 
Zel'dovich (tSZ) and kinetic Sunyaev-Zel'dovich (kSZ) effect statistics. The results derived here are complementary to calculations based on the 
halo model and are applicable to tSZ effect originating from the extragalactic ionized medium. 



3 MIXED LOWER ORDER MOMENTS: CUMULANTS AND CUMULANT CORRELATORS 

In this section, we present results in real space and in the harmonic domain that are completely generic and can provide useful information to 
study the pressure fluctuations in the baryonic gas. We apply these results to understand the studies involving tomographic bins from weak lensing 
surveys. In case of CMB lensing one of the source planes is identified with the last scattering surface. The cross-spectra P^^ (defined below) can 
be probed by cross-correlating tomographic weak lensing maps kj;) and the tSZ y maps of projecte d pressure flu ctuations. These correlations 
sample the three dimensional density-pressure cross-spectra Pa^ ; using the small angle approximation ( lKaisej[l998l) one can write: 

,r hi{es)WTH{ies); P^s{k,z) =b^{z)Ps{k,z). (8) 

Here bi(6s) and Wtu{103) ~ (2Ji {I6s)/I9s) are Gaussian and Top-hat windows respectively. The line-of-sight directions Oi and Q.2 are separated 
by an angle 6^12. In our notation |1| = /. The weight function ajf^j (r) for tomographic convergence and ujsz{r) for tSZ surveys are defined in Eq. l|7j 
and Eq. l|4j respectively. 

To study non-Gaussianity in pressure fluctuation we need to go beyond the power spectrum analysis; we propose the mixed cumulant corre- 
lators for this purpose. Cumulant correlators have the advantage of being very simple to estimate, and are defined in real space so can be useful 

© 0000 RAS, MNRAS 000, 000-000 
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Figure 4. The projected (or 2D) weak lensing and tSZ cross-spectra is plotted as a function of harmonics I (left panel). We also show the projected weak lensing 
convergence power spectra as a function of I (right panel). Two different values of erg = [0.8, 1.0] are considered. 



for smaller surveys. Similar results can be obtained for the kurt-spectra, which we will not consider here. We will consider a top-hat smoothing 
WTiiilds) for the convergence maps and a Gaussian beam bt{9s) for the 1/3(^1) maps. 



(10) 



Here B^-^^ represents the mixed bispectrum involving three dimensional density contrast 5(x) and pressure fluctuation 7r(x). These mixed cu- 
mulant correlators should be considered in addition to the pure ones that probe and B^^^. These results can trivially be extended to higher 
order to compute mixed cumulant correlators. At the level of fourth order they will probe the mixed trispectra of various types; e.g. we can use 
{K^i)(ni)j/s (f22))c (K(i)(f2i)2/s(02))c to probe the trispectra T^^^^, T^^^^ and Tsir-Kir- At these level these statistics are completely general 
and they can be estimated in real space simply by cross-coorelating smoothed y{9s) and k{9s) maps raised to suitable powers. The higher order 
expressions follow from the same logic. 

The derivation so far has been completely generic . It does not depend on specific form of the multispectra. Later we will be use a specific form 
for the correlation hierarchy which will allow us to include correlation functions to all order. This is done using a generating function formalism. 

The results in the real space are suitable for smaller surveys. However with the increase in survey size it is pertinent to consider a har- 
monic space approach as real space measurements are often correlated. We will introduce power spectra in the harmonic domain that are Fourier 
transforms of the corresponding cumulant correlators in real space. 



rnn.y / (2?! + 1) (2^2 + 1) / h h I \j^..y. rW-^^X^ / (2Zl + 1) (2/2 + 1 ) ( h h I \ ^yy. 

' 47r(2/ + l) V "^2 / ' 47r(2Z + l) V "^2 m ) 'i'^'' 

The skew- and kur t-spectra are simple to estimate from real data and the scatter for such estimates is well understood l lMunshi & Heavensll201(il ; 
iMunshi et alj201 0) . These results can be used to construct a family of skew-spectra for various choices of tomographic slices. We also provide the 
expressions for the auto- and cross-spectra for tSZ and wl surveys below: 

Jo d\(r) \dA{r) > 

The auto spectra for wl and tSZ can be expressed in terms of P^^ and V-,,-^ with suitable changes in the weight functions, i.e. (r) for wl surveys 
and iJl^ (r) for tSZ surveys. 



3.1 Hierarchical Ansatze 

In deriving the above expressions we have not used any specific form for the matter correlation hierarchy, however the length scales involved in 
small angles are in the highly non-linear regime. Assuming a tree model for the matter correlation hierarchy in the highly non-linear regime one 
can write the most general case as (White 1979; Peebles 1983; Fry 1984; Bernardeau & Schaeffer 1992; Szapudi & Szalay 1993; Bemardeau & 
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Figure 5. The bias function b{r]) of the reduced convergence is depicted for two different smoothing angular scales. The dashed lines show the analytical predictions 
from the hierarchical ansatz and the solid lines show predictions from the lognormal distribution. The left panel corresponds to 6s = 10' and the right panel 
corresponds to 8s = 5'. The different curves in each panel correspond to different redshift bins. The lower curves correspond to shallower (lower redshift) bins and 
the higher curves correspond to deeper (higher redshift) bins. The top-most curve coiTesponds to the projected survey without any tomographic information. 




Schaeffer 1999): 

(JV-l) 

^^(ri,...rN) = Q'^'" n ^2('"'''"j)- (13) 

a, N — trees labellings edges(i,j) 

It is interesting to note that an exactly similar hierarchy develops in the quasi-linear regime in the limit of vanishing variance (Bernardeau 1992), 
however the hierarchal amplitudes QN,a become shape-dependent in such a case. In the highly nonlinear regime there are some indications that 
these functions become independent of shape parameters as has been suggested by studies of the lowest order parameter Qa = Q using high 
resolution numerical simulations (Sccociamarro et al. 1998). In the Fourier space such an ansatz will mean that the whole hierarchy of multi- 
spectra B^, Tg can be written in terms of sum of products of power-spectra Pg, e.g. in low orders we can write: 

Bi(ki,k2,k3)2fc.=0 = Q{Ps{ki)Ps{k2) + Ps{k2)Ps{h)+ Ps{k3)Ps{ki))5D0il +1^2 + 1^3), (14) 

Ti(ki,k2,k3,k4)5;fe^=o = [iiaPi(fci)Pi(lki +k2|)Pi(|ki +k2 +k3l) +cyc.perm. 

+RbPs{ki)Ps{k2)Ps{k3) + cyc.perm.](5i3(ki + kz + kg + k^) (15) 

Different hierarchal models differ in the way they predict the amplitudes of different tree topologies. Bernardeau & Schaeffer (1992) considered the 
case where amplitudes in general are factorizable, at each order one has a new "star" amplitude and higher order "snake" and "hybrid" amplitudes 
are constructed from lower order "star" amplitudes (see Munshi, Melott & Coles 1999a,b,c for a detailed description). In models proposed by 
Szapudi & Szalay (1993) it is assumed that all hierarchal amplitudes of a given order are actually degenerate. We do not use any of these specific 
models for clustering and only assume the hierarchical nature of the higher order correlation functions. In the past, primarily data from galaxy 
surveys have been analyzed extensively using these ansatze. Our main motivation here is to show that cross-correlation statistics of weak-lensing 
surveys and the tSZ surveys can also be analyzed using such techniques. The most general form for the lower order cumulant correlators in the 



large separation limit can be expressed as: 

{til{(li)5ys{(l2))c = 2Q3C3[Xe,2:e,,] =C'27'C3[2e,XeiJ = C'27(fc^)e(5y.(ni)fc,(n2))c, (16) 

{Kt{h^)5ys{Q.2))c = (3i?a+6i?6)C4[Xe,/9iJ =C';?i'''C4[X|3/£.iJ =C37{^?)^{K4J7i)5y4n2))c, (17) 

{Kl{Cl^)5yl{Sl2))c = mbCi[Tlle,^]=Cll^'c,[Ie^Ie,,] = C^ni^l)c{5yl)c{iis{^i)^^ (18) 

{4{Q.i)5ys{(l2))o = (245. + 36S5 + 45^)4 Pl.IeiJ = QT'CsPI^X^i,] = C^^ &{i^s{(li)5ys{(l2))o, (19) 

{^s{(li)5yl{(l2))c = (125„ + 6S6)C5[/|,/fliJ = CV,:^' Cs[lll0rA = C;^{i4)l{5yl)o{^s{(li)5ys{(l2))c. (20) 

Here Cpq denotes the cumulant correlators of the convergence field and Cpq denotes the cumulant correlators for the underlying mass distribution. 

The subscript c denotes the connected components of the higher order correlation functions. The amplitudes Ra — 1^2, Rb ~ I's and Sa = 



vij Sb = v:iV2-, Sc = 1^4 are expressed in terms of the vertices which can be evaluated using HEPT (Hyper Extended Perturbation Theory; 
IScoccimarro et al. I( ll998h ). The exact expressions for various Cpg presented in Eq. ilOl are obtained by counting of the relevant tree diagrams. In 
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PDF and Bias of Reduced y Parameter 

1 10 h H 




1 2 3 1 2 

(5y 6y 



Figure 6. The PDF (left panel) and bias function (right panel) b(> 5y/ (y)) are plotted as a function of Sy/ (y). They are computed using two different approaches. 
The dashed lines show the analytical predictions from the hierarchical ansatz and the solid lines show the predictions from the lognormal distribution. The lines 
correspond to different angular scales as depicted. 



the limiting situation when fli = fl^ we can recover the corresponding cumulants. Extending the above results to arbitrary order we can write: 

This is a generalization of the usual definition of cumulant correlators for the case of two different fields in this particular case convergence k 
and Sys- The smoothing angular scales and the window function are left completely arbitrary. These definitions can be also used to define mixed 
normalised one-point S'^ parameters (with N — p + q) involving two different fields by considering the limiting situation Cli = Cl2. 

c;+, [[lo.nr.] - '"-'P,,f,^l^,^^^\ P^Jw, = / (^) wUWs); (22) 

Pdsz ^ / (^) . / (-^) W..mHO'.)eM^-0..). (23) 

Though results can be derived for different smoothing angular scales, for simplicity we will only consider identical smoothing beam size 6s = Og. 
The hierarchical expression for the lowest order cumulant i.e. S3 for convergence was derived by Hui (1998). He showed that his result agrees well 
with numerical ray tracing experiments. Later studies have shown that higher order cumulants and even the two-point statistics suc h as cumulant 
correla tors can also be reliably modelled in a similar way (Munshi & Coles 1999; Munshi & Jain 1999a). More recently it was shown jMunshi et al.l 
l201lj) that the statistics of tSZ too can be modelled according to the same prescription. In particular it was shown that the lognormal distribution 
can be used to predict the PDF and bias associated with the tSZ maps. We extend these results to the case of joint analysis of weak lensing and tSZ 
maps. 

We will develop these results further to construct the full joint 2PDF simply using the individual bias functions for tSZ and weak lensing. The 
hierarchical ansatz allows us to write the joint 2PDF as: 

and its relation to the bias associated with collapsed objects in underlying density field rj = 1 + S. 

As an aside, it is simple to check that if we are dealing with tomographic redshift slices and Ks''' of the same or different weak lensing 
surveys the corresponding cumulant correlators are defined as: 

r^s iJ" (r)Lo'' (rW(r) 

C;^, [[2:«J-^[X.jr^[X«.J] = / ''' Jj:l\'/, ' [IsAiV'[ToAir''[Te,.] dr. (25) 

^0 (r) 

The expressions for the cumulant correlators are constructed by replacing ys{^) in Eq. ( I21t with K,i''\Cl): 

(4,(^1)4,(02)), = c.Vp+^fPej'^+'-'^p.iJl = c;,^(4))(''-''(4))i'-''(4''(^^i)4^'(j^2))c. (26) 

These expressions are quoted here for the highly nonlinear regime. As is well known a similar hierarchy developes in the quasilinear regime 
so the same analytical tools are applicable in such a situation. We will next use these expressions to compute the lower order moment of one and 
two-point PDFs. 

The corresponding results for the Ostriker-Vishniac effect (and kinetic Sunyaev-Zel'dovich (kSZ) effect; ICastrd ( l2003h ) will be presented 
elsewhere. 
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4 JOINT PROBABILITY DISTRIBUTION FUNCTION FOR k.s AND ys 

To compute the bias associated with peaks in the convergence field we have to first develop an analytic expression for the generating function 
l3Ky{zi, 22) for the convergence field Hs and the tSZ field Sys = ys ~ iVs}- For that we will use the usual definition for the mixed two-point 
cumulant correlator Cpg : 

{fci)r'{5yl)r'{^.(ni)<5y.(f^2))c ' 
(for a more detailed description of cumulant collators see Munshi & Coles, 1999b). We will show that, like its density field counterpart, the 
two-point generating function for the convergence field Ks can also be expressed (under certain simplifying assumptions) as a product of two 
generating functions I3{z) which can then be directly related to the bias associated with "hot-spots" in the convergence field. It is clear that the 
factorization of generating function actually depends on the factorization property of the cumulant correlators i.e. Cpg = Cp^Cgi . Note that such 

a factorization is possible when the correlation of two patches in the directions Qi and O2 {his{Cli)5ys(Cl2))c is smaller compared to the variance 
for the smoothed patches: 

Here zi and Z2 are dummy variables. We will now use the integral expression for cumulant correlators (Munshi & Coles 1999a) to express the 
generating function which in turn uses the hierarchical ansatz and the far field approximation as explained above: 

It is possible to further simplify the above expression by separating the summation over dummy variables 21 and 22, which will be useful to 
establish the factorization property of two-point generating function for bias fi^yiz-i, 22). We can now decompose the double sum over the two 
indices into two separate sums over individual indices. We do not use any of these specific models for clustering and only assume the hierarchal 
nature of the higher order correlation functions. In the past, primarily data from galaxy surveys have been analysed extensively using these ansatze. 
The motivation here is to show that cross-correlaion statistics of weak-lensing surveys against the tSZ surveys can also be analysed using such 
techniques. The most general result for the lower order cumulant correlators in the large separation limit can be expressed as: 

R - r^rlr-^^ (r-\^'^^2 S!^^fl / ^^^wlW FT- , \ i^Mfh R ( ^Sz(r-)b^(r) \ 

/3.„(2„22)-y^ drd,W^^/?.(^^[X«Jwtj|^/?v(^^r(;:^P«Je^j- (30) 

The above expression is quite general and depends only on the small angle approximation l lLimbeiiri954l) and the large separation approximation 
and is valid for any given specific model for the generating function Q^ij). However, it is easy to notice that the projection effects as encoded 
in the line of sight integration do not allow us to write down the two-point generating function jS^yizi, 22) simply as a product of two one-point 
generating functions /3k (21) and I3y{z2), as generally is the case for the density field r] — 1 + 5. 

As in the case of the derivation of the probability distribution function for the smoothed convergence field Hs , it will be much easier if we 
define a reduced smoothed convergence field ria . The statistical properties of rjs are very similar to that of the underlying 3D density field (under 
certain simplifying approximation) and are roughly independent of the background geometry and dynamics of the universe. In a similar manner we 
also define the reduced ys field 77^. We define the reduced convergence ris and the reduced tSZ y map rj's respectively by the following expressions: 
rjs = (ks — «;m)/(— Km) = (1 + Ks/\nm\) and rj's — Sys/{ys)- Where the minimum vale of Km is defined as: km = — J^" dr u)^\(r-) and 
{y) is the average value of the y parameter. It is easy to notice that the minimum value of the convergence field will occur in those lines of sight 
which are completely devoid of any matter, i.e. 5 = —1 all along the line of sight. We will also find out later that the cosmological dependence 
of the statistics of the Hs field is encoded in km and this choice of the new variable 77s makes its related statistics almost independent of the 
background cosmology. Repeating the above analysis again for the 77s field, we can express the cumulant correlator generating function for the 
reduced convergence field 77s as: 

While the above expression is indeed very accurate and relates the generating function of the density field with that of the convergence field, it 
is difficult to handle for any realistic practical applications. Also it is important to notice that the scaling functions such as h{x) for the density 
probability distribution function and h{x) for the bias associated with over-dense objects (with a:: = (1 + 5)/ ^2) typically estimated from 
numerical simulations especially in the highly non-linear regime (^2 is the volume average of two-point correlation function). Such estimations 
are plagued by several uncertainties, such as the finite size of the simulation box. It was noted in earlier studies that such uncertainties lead to only 
a rather approximate estimation of h{x). The estimation of the scaling function associated with the bias i.e. b[x) is even more complicated due to 
the fact that the two-point quantities such as the cumulant collators and the bias are more affected by the finite size of the catalogs. So it is not 
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fruitful to actually integrate the exact integral expression we have derived above and we will replace all line of sight integrals with its approximate 
values. The previous study by Munshi & Jain (1999) have used an exactly similar approximation to simplify the one-point probability distribution 
function for Ks and found good agreement with ray tracing simulations. We will show that our approximation reproduces the numerical results 
quite accurately for a wide range of smoothing angle, 

\i^m\ ~ ]^rsU]„i{rc)\ (y) ~ irsa;sz(rc)6,r(rc); < re < r^; (32) 
/ 2\ _ 1 ^li{rc) \ f dH I . 2 f,a \\ /X 2\ _ 1 bl{rc)iOsz{rc) \ f (fl , I . 2,0 J ,0,-, 

{n^{Cl,)Sy^{Cl2)}. ^ r° '^''d>r)'^''^^ [ / ^P^(j^))WTu(iesMe.)e^p[i\-e^2]\. (34) 

Use of these approximations gives us the leading order contributions to these integrals and we can check that to this order we recover the factor- 
ization property of the generating function i.e. (zi, 22) = /?,j(2i)/3^/(z2), 

/3,,,,'(2l, 22) = l3r,{zi)P^,{z2) = Pl + 5{zi)l3i + S'{z2) = t{zi)t{z2). (35) 

So it is clear that at this level of approximation, due to the factorization property of the cumulant correlators, the bias function brj{x) associated 
with the peaks in the convergence field Ks , beyond certain threshold, obeys a similar factorization property too, which is exactly the same as its 
density field counterpart. Earlier studies have established such a correspondence between the convergence field and the density field in the case of 
the one-point probability distribution function p{S) (Munshi & Jain 1999b), 

bj^(xi)hj^(xi)b^, (x2)h^> (x2) = bi+sixi)hi+sixi)bi+g, (x2)hi+s, (x2)- (36) 

The following relation between Pn{z) and 6^ (a;) is useful for modelling. However, for all practical purpose we find that the differential bias /3>,, (y ) 
as defined is simpler to measure from numerical simulations due to its cumulative nature: 

1 r°° 1 r°° r(z) 

br,{x)hn{x) — / dzr(2) exp(a;2); bn{> x)hn{> x) — / dz exp(3:2). (37) 

Notice that the fuction t{z) defined as t{z) — z^Cpi/p\ also acts as a generator for the Cpi parameters. The entire hierarchy of the Cpq can 
be constructed from the Cpi parameters i.e. Cpq — CpiCqi at large separation limit. Using these results we can finally write down: 

Pvv' iv, v')dvdv' ^Pviv)Pv' + bv{v)^i2 bn'{v'))dvdv'; V ^ ks/\km\; r)' = Sys/ {ys); ^"^^ = {ri{{li)r]' {CI2)) . (38) 

It is important to notice that although the bias brj(x) associated with the reduced convergence field and the underlying density field are exactly 
same, the variance associated with the density field is very high but the projection effects in the convergence field bring down the variance in the 
convergence field to less than unity which indicates that we have to use the integral definition of bias to recover it from its generating function (see 
Appendix A for more details). Finally, writing down the joint probability distribution function p^y {Ks,ys) for the smoothed projected convergence 
field Us and the ys in terms of their reduced versions ri and 77' : 

Pi^y{Ks,ys)dKdy = pn{K,s)py{ys){l + &k(ks) Cil by{ys))dKs dys\ 

hy{ys) = b^(ys/{ys))/{ys); 6«(ac,) = 6^(^c,/|^,„|)/|Ac,„i; C^l = {^^{(1^(12)) = iTi' KMys))- (39) 
In an earlier study, Munshi et al. (1999b) used similar agruments for the convergence maps to show thatpK(Ks) = p-i^[k.s /\i<.m\) /\km \ and recently 
it was generalised by (Munshi et al. 20 1 1 ) for the y-maps to py (ys ) ~ p-qiys / {ys)) / (ys) if we notice that — I ( I '^m I (2/3 ) ) ; then the above 
expressions helps us to write b^(Ks) ~ bn{Ks/\Km\)/\nm. \ and by{ys) — b^{ys/ (ys))/ (ys) This is one of the main result of our analysis. 

A similar result can of course be derived for cross-correlation of tSZ ys map and weak lensing convergence maps k^'^ from individual 
tomographic bins. The joint PDF p'-'' {ni'^ , ys) for tomographic maps /ti*' and projected ys maps can be expressed in terms of the individual PDF 
p^^\K.s),p{ys) maps and the bias &re'(K:s) and by{ys): 

p^'\K[^\ys)dnf dys =p«(.c.)p(y.)(l + &«(4'') 

= fe,(r,)/|fc«l; by{ys) = b^{r^') / {ys); = e^a'Vd'^^' K^.))- (40) 

Notice that the bias and one point PDF for both the projected convergence maps k,s and tomographic maps t^k^ as well as for the tSZ map y are all 
constructed from the same underlying PDF Pn and the bias 6,, associated with the 3D density distribution. The individual PDF Pk.(i^s) and py (ys ) 
and the bias b^ii^s) and by(ys) have already been tested against numerical simulations and were found to be remarkably successful. 

This technique will also allow for the computing of joint PDF of convergence maps for CMB lensing and the tSZ ys maps. This can be 
achieved by replacing the k„i in Eq. J39| l with corresponding him for the last scattering surface (LSS), i.e. tim = — j^^^^ drw^\{r) where the 
source distance Ts is now replaced b y the comoving distance to the LS S tlss- This is especially relevant given the recent evidence for a CMB 
lensing signal with WMAP and ACT jDas et aLll201 iklSmidt et al.ll201 ih The results for the cumulant correlators can be modified analogously. 
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Though the resuUs Eq. l |39l ) and Eq. ( I40t are derived using the hierachical ansatz the final results are remarkably independent of details of 
any of the assumptions that were used in deriving them. This is an indication of their more genric validity. Indeed it has been shown that the PDF 
and the bias for the tSZ and weak lensing fields can be constructed from equivalent but other valid descriptions of PDF and bias of underlying 3D 
density constrast. The lognormal distribution for the underlying PDF and bias provide one such description. In Appendix-B we have provided a 
short description of the lognormal distribution. The lognormal model based PDF and bias for and ys are given by: 

Pk(ks) = Pin(Ks/|K,n!)/|fc„|; Py{5y s) = Pin{5y s / {y s)) / {y s) \ = fein(K/|«:m|)/|fc„,l; by{5y s) ^ bi^{5y s / {y s)) / {v a) ■ (41) 

The resulting PDF and bias matches with the ones from perturbative calculations for higher smoothing angular scales. 

As mentioned earlier, the joint probability of pKy (acs ,ys) is a noisy statistics. The integrated measure or cumulative PDF of the convergence 
crossing a threshold Hs and the tSZ field ys crossing a threshold kq and yo respectively is given by: 



PKy{> fio,> yo) = / dhis / dy p^y{tzs,ys); Pk{> ^o) 



dKsPK{Ks); Py{> yo) 



dyspijjs); 



/oo poo poo roc 

dtisP[^s)hti{His) I I dtisp{i^s); by{> yo) = / dhisPy[ys)hy{ys) / I dysPy{ys)- 



(42) 
(43) 



The bias b{> kq) is plotted in figure (O for different redshift bins. As is clear the high k, regions in deeper redshift bins are more biased compared 
to the lower redshift bins. The bias is higher for larger smoothing angular scales. The figure ([SJl shows corresponding plots for the tSZ surveys. 
The following statistics can be constructed from the cumulative bias to probe the joint SZ and weak-lensing cross-correlation to all orders: 



B{> K.0, > yo) = [6k(> K.o)by{> yo)]^^'^ = 



PKy(> KQ, > yo) 



Pk(> l^o)Pyi> yo) 



1/2 



(44) 



It is important to notice that our modelling of the PDF p^ and py is independent of the modelling of respective variances. The variance calculations 
for ys depend on inputs such as the detailed modelling of bias and its redshift evolution that can be independently checked by comparing with 
simulations. However, given a correct variance the lognormal model or the hierarchical ansatz can be used to model the entire PDF as we have 
shown. For construction of the PDF of the scaled or reduced tSZ i.e. 77 we also need the average (ys)- The modelling of average (ys) is independent 
of the construction of PDF. However what we have shown here is that given these two inputs we can reliably predict the PDF of y. While the PDF 
of reduced tSZ field rj is independent of cosmology the dependence of y parameter is encoded in the definition of reduced ys field 77' = Sys/ (j/s). 
The convergence PDF is defined in terms of the PDF for the reduced convergence 77 — «:,„]. The PDF for the reduced convergence Ks as well 
as the scaled tSZ map ys are both shown to be same as the PDF of underlying ark matter distribution 5. The dependence of on cosmological 
parameters is encoded in the definition of K,n. The cosmological parameter dependence of bias b^i^s) and by{ya) is encoded in (ys) and Km- 

It is also interesting that our calculations show that the PDF of projected or 2D con vergence as well as the tSZ maps are described by the 
underlying 3D density contrast 5 which is in agreement with l lMunshi & Jaii]l2000l[200lh and provides robust mathematical justification for the 
use of lognormal model and its link to linear biasing prescription. 

A simple order of magnitude estimate for the weak lensing cumulants is Sp = )/ {ti^s)^^^ is given by — Sp/\KmY'^'^ ■ The cumulants 
of the field Sys/ {ys) simply follows that of the underlying density constrast. The normalised cumulant correlators of the smoothed convergence 



field K3 and the field Sys/ {ys) is given by Cp| 



ip-i 



5 INCLUSION OF NOISE 

In our discussion so far, we have ignored the effect of noise in the estimation of one- and two-point PDFs. If we assume the noise n to be 
uncorrrelated with the signal, the noisy one-point PDFs can be expressed as a convolution of the signal and noise PDF: 

/oc poo 
pl^lif^s -n)gi'^{n)dn; p„iys) = py{ys - n)gy{n)dn. (45) 
- oc J —00 

We have denoted the i-th tomographic noisy convergence PDF by The noisy J/-PDF is denoted by p„iy)- The noise in y-map is due the 

detector noise and in k maps due to intrinsic ellipticity distribution of galaxies. The noise PDF in each case is often assumed to be Gaussian; we 
will use (/k ' (n) and Qy (n) to denote the PDFs of the noise in y-map and tomographic k maps respectively. The modification of the joint PDF of 
Ks and ys due to noise can be computed using a similar convolution: 



/OO 
P^''\ks — ni,ys — 712) g'i/\ni) 32,(712) d7ii dn2 
- 00 



(46) 



Using Eq. ((39) with Eq. l |46l l, we derive the following expression: 

n'i\Ks,ys)=P^^{Ks)pM{^ + b'i\Ks) h^iVs))- (47) 
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The inclusion of noise tlius preserves tlie factorization of thie 2PDF whicli can be expressed in terms of noisy PDFs and bias functions. The presence 
of noise does not affect the two-point correlation function ^J'J' ; The noisy bias functions b^^ (/t) and for k and y are respectively defined by 

the following expressions; 

/oo /'OO 
- n)&W(Ac, - n)g^{n)dn / / p«(ft, - n)g'i\n)dn (48) 
-ioo J — oo 

/oo POO 
p'i\Ks)h'i\Ks)(n)sdKs / p'i\t,)dK- (49) 



We have defined a cumulative bias function fen' (> k) for the noisy map. It is constructed from the bias fel'' (k) and the one-point PDF pi*' (k) in an 
exactly similar way as the noisfree PDF and bias function Eq. {43}. The noise in convergence maps depends on the intrinsic ellipticity distribution 
of galaxies o^. For a smoothing radius 9s the noise variance Ok that completely determines the noise PDF gi' '(ni) is given by: Gk. = /2img9^. 
Here Ug is the number density of galaxies per arcmin'^. Similar expressions can be used to define the fe^(> y): 



KiVs) = I Pyivs - n)by{ys - n)gyjn)dn / pl^^ {ks - n)gl''' {n)dn (50) 

J —ioo J — oo 

poo poo 

bn{> ys) = p'i\ys)b\l\ys)dys / p„(ys)dys. (51) 

These expressions are sufficient to compute the cumulative bias for weak lensing tomographic slices and y-maps; these bias functions are enough 
to construct the entire family of two-point cumulant correlators involving weak lensing and tSZ maps even in the presence of noise. 



6 CONCLUSION 

We have presented a detailed analysis of the higher-order (non-Gaussian) cross-coiTelation of tSZ maps and the projected as well as redshift- 
resolved tomographic maps from weak lensing surveys. 

What distinguishes our work in this paper is that it extends earlier studies which have focused mainly on cross-correlation analysis at the 
level of power spectrum. The tSZ maps as well as the maps derived from weak lensing surveys are intrinsically non-Gaussian even at large angular 
scales. In this study we extend the cross-correlation statistics to higher-order. We employ the mixed cumulants and their correlators that can extract 
higher-order statistics that are useful in probing the thermal history of the Universe when studied in association with redshift resolved tomographic 
weak lensing data. The statistics that we introduce can be used in projection (2D) or with tomographic information, independent of any analytical 
modelling to analyze simulations or observational data. The estimators presented here are generic and can be used for other cross-correlation 
studies involving two different arbitrary fields. The mixed cumulants and their correlators defined in this paper are the lower order moment s of 
the corresponding one-point or two-point PDFs an d associated bias. These statistics generalize cumulant correlators defined in l lMunshillioOol) for 
weak lensing surveys and in ( lMunshietal]|20Ilallbh for tSZ datasets, and are useful in cross-correlating these two independent data sets. 

Beyond the order-by-order calculation, a generic prediction for the cumulants or their correlators to an arbitrary order is required. This can 
be achieved by adopting the generating function formalism inherent in models based on the hierarchical ansatz. We have employed a generic 
version of the hierarchical ansatz as well as the lognormal distribution to model the underlying mass distribution. Both of these have been studied 
extensively in the literature. These particular models go beyond the lowest order in correlation hierarchy. For the hierarchical model, the generating 
function is parametrized in terms of a single parameter uj which is fixed using inputs from numerical simulations, and is the only freedom allowed 
in this model. In the perturbative regime this parameter can be calculated analytically. Indeed, modelling of pressure fluctuations in terms of the 
underlying den sity distribution is mo re complicated, and we rely on a redshift dependent linear biasing model that is expected to be valid at large 
angular scales jMunshi et all 20 II al ibi). More complicated models can be considered in this framework using additional inputs from numerical 
simulations; here we have tried to keep the modelling as simple as possible. The halo model can be used for the construction of the lowest order 
cumulant correlators. However, in such an approach, results can be derived only in an order-by-order manner and the relevant PDF and bias needs 
to be constructed using suitably truncated Edgeworth expansions. As mentioned previously, the PDF and bias for the the weak lensing data have 
been studied by many authors using both the hierarchical ansatz and the lognormal distribution. The studies have shown how accurately these 
distributions can model numerical simulations. In more recent studies, the hierarchical ansatz has also been used successfully in tSZ modelling. 
Based on this success, our primary goal in this study has been to extend such analysis to compute the joint PDF of and j/s along same or 
different lines of sight. In lieu of individual haloes, we consider here the diffuse tSZ component that correlates with the large-scale weak lensing 
shear and convergence. 

We have shown that the two different variables, i.e. the scaled tSZ parameter 77' = Sya/{ys) and the reduced convergence field 
can be introduced to simplify the analysis. Employing the hierarchical ansatz, we next show that the statistics of these reduced variables are the 
same as the underlying density p/pb = 1 + 5 under certain simplifying approximations, which is true both at the level of one and two-point 
statistics. Thus, the final results do not depend on the specific details of the input parameters of hierarchical ansatz, and may point to a wider 
applicability. The joint distribution involving Ks and ys can thereby be written in terms of p{i]) and b{rj) with suitable scaling involving (ys) and 
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Hm ■ The detailed construction of the PDF is insensitive to modelling of the bias of the hot ionized gas with respect to the dark matter distribution. 
The detailed prescription for such bias enters the results through the computation of the variance and the two-point cor relation function ^iik(^12)- 

Specific models of hierarchical clustering have prev iously been employed to understand the tSZ effect (see e.g. dValageas. Schaefer & SiM 
l200ll ; l\^lageas & Silldll999l : IValageas & Schaeffeill2000h ). However, they have been employed to model contribution from the collapsed objects 
(e.g. massive clusters). In this case, the tail of the scaling functions h{x) and b[x) are relevant as they describe the collapsed objects. In addition 
to these inputs, the modelling would involve a prescription for the hydrodynamic equilibrium of the gas residing in the halo. These can be used 
to model the PDF of yip,), i.e. Py(y) or by{y). The modelling considered here is complementary to our analysis. We have considered redshift 
dependent linear biasing to model the shock heated gas in the inter-galactic medium that produces a diffuse tSZ effect. The linear biasing model 
has been studied at the level of the power spe ctrum by many author s (Goldberg & Spergel 1999a b; Peiris & Spergel 200 0j). The hier archical model 
that we use here was also previously used bv lZhang. Pen & Wand J2002f) . I Coorav. Tegmark & Hul ( 12000) and ( ICooravll2000l. l200l') for modelling 
lower order moments where the linear biasing model for the pressure fluctuations was assumed. Our results are natural generalizations of such 
studies to higher order multi-spectra to study the non-Gaussian aspects of tSZ maps. 

The generic statistics that we have developed for lower order cumulant and their correlators can be used to probe the higher order correlation 
functions in a model independent way. The redshift-resolved tomographic weak lensing maps when cross-correlated with tSZ maps can provide 
valuable information of the redshift evolution of the clustering of the ionized gas. The cumulant correlators are collapsed higher order correlation 
functions and typically carry more signal to noise than ordinary higher order correlation functions. Hence, they are useful for analyzing data from 
numerical simulations and observational maps. We have presented detailed analytical derivations regarding the inclusion of noise in our analysis 
that follow from simple convolution. It is important to realize that the modelling we have presented is highly constrained without much room for 
any additional fine-tuning. The bias that appears in the independent analysis of the tSZ maps (i.e. by{ys)) and the bias that is required to model the 
statistics of weak lensing maps (i.e. (ks)) are enough to model the higher order cross-correlations between these two maps to an arbitrary order. 
This could provide a much needed consistency check for any possible cross-contamination, as well as inputs from "gastrophysics". It is moreover 
interesting that the tail end of the same scaling functions h{x) and b{x) that we have adopted here are used to model the tSZ. 

It is known that the IGM is most likely have been preheated by non-gravitational sources. The feedback fro m SN or AGN c an play an impor tant 
role. The analytical modelling of such non-gravitational proce sses is rather difficult. Numerical simulations dSpringel etail l2001: da Silva e t al.l 
l2000ll2004lwhite, Hernguist & Springell20^ : lLin et alj2004l) have shown that the amplitude of the tSZ signal is sensitive to the non-gravitational 
processes, e.g. the amount of radiative cooling and energy feedback. It is also not straightforward to disentangle contributions from competing 
processes. The inputs from simulations are vital for any progress. Our analytical results should be treated as a first step in this direction. We have 
focused mainly on large angular scales where we expect the gravitational process to dominate and such effects to be min imal. Indeed, stud ies 
involving simulations have found that the cross-spectra involving tSZ and wl maps are less affected by the baryonic physics dShao et al]|201 lal) . 

Our results are derived using cross-correlating tomographic slices from weak lensing surveys and the projected tSZ surveys. With suitable 
modification, they can be equally applicable to cross-correlation studies using tomographic maps from the same or different weak lensing surveys. 
The results can also be used to study the weak lensing of CMB by cross-correlating it against tomographic weak lensing maps or tSZ maps. To 
move beyond a tomographic or projected s urvey, a complete 3D analysis can be invoked. An analysis similar to what has been presented here for 
the kinetic Sunyaev-Zel'dovich effect (e.g. Ishao et ail hoilbl) ) can provide valuable information about the reionization history of the Universe 
(Munshi et al. 201 1; in preparation). 
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APPENDIX A: GRAVITATIONAL CLUSTERING AND THE HIERARCHICAL ANSATZ 

The lower order cumulants and cu mulant correlators as well as the one- and two-pint PDFs are commonly used statistics in cosmology to qunatify 
clustering dSernardeau et all 1 20021). Use of generating function techniques t o go beyond order by order calculations is an im portant milestone 
in this direction jPeeblesI ll97lUBalian & Schaeffej|l989l : lBernardeai]|l992l : lBemardeau & SchaeffeJI 19921: iBernardeatjl 19941) . This approach is 
complimentary to the halo model based approach that rely on Press-Sechter mass formalism i Press & Schechtej 1974 ). We provide a very brief 
review of this formalism in this appendix. 



Al The Generating Function 

In scaling analysis of the probability distribution function (PDF) the void probability distribution functi on (VPF), denoted as P,,(0), plays mo st 
fundamental role, which can related to the generating function of the cumulants or Sn parameters, 4>{z) ( IWhitell979l ; lBalian & Schaeffe jl 19891) : 

P40) = exp(-iVa(iV,)) = exp ( - <?!>(iVc)/Cl) • (Al) 

Where Pv(Q) is the probability of having no "particles" in a cell of of volume v, N is the average occupancy of these "particles" and Nc ~ iV^I 
and (j{z) = —(j>{z)/z. Here (j>{z) = X^^i Sp/p\ is a generating function for ^jv parameters. The genrating function Gjr) for th e vertex 
amplitudes Up is related to 4>{z)- A more specific model for G{t) is generaly used to make more specific predictions jBalian & Schaeffe3l 19891) : 

Gs{t)= (^l + r/n^y'"'. (A2) 

The parameter Kg is related to the parameter uja to be introduced later. The two generating function are related by the following expression 
tealian & Schaeffer.il989, : ,Bernardeau & Schaeffer 1992, 1999j): 

0(2) = zQs{t) — ^zt-^C;s{t); t = —z-^Qs{t). (A3) 
/ dr dr 

However a more detailed analysis is needed to include the effect of corr elation between two or more correlated volume element whi ch 
will provide information abou t bias, cumulants and cumulant correlators l lBemardeau & SchaeffeJ 1 19921 : iMunshi. Coles & Melotj Il999allbt 
iMunshi. Melott & Coles! [T999ch Notice that t{z) (also denoted by I5[z) in the literature) plays the role of generating function for factorized 
cumulant correlators C;^' (Cp, = C^^^Cy ): t{z) = E^^i C^iVp! z^ . 

A2 Generating Functions and the Construction of the PDF and Bias 

The hierarchical form of higher order correlation functions appear in two completely different regimes in gravitational clustering. The generating 
function approach is conveninent technique to sum up to arbitrary information that leads to construction of entire probability distribution functions. 

A2. 1 The Highly Non-linear Regime 

The PDF p{S) and bias h(5) can be related t o their generating functions VPF 4>{z) w\At{z) respectively by following equations jBalian & Schaeffeil 
ll989l : lBemardeau & Schaeffe3ll992l[l999l) 

r°° dz \{l + S)z-(|,{z)^ ,,,,,,, r°° dz r(l + 5)2-0(2)] 

It is clear that the function (j){z) completely determines the behavior of the PDFp((5) for all values of S. However different asymptotic expressions 
of (j>{z) govern the behavior of p{5) for different intervals of 5. For large y we can express (f){z) as: (j){z) = az^^". Here we have introduced a new 
parameter u for the description of VPF. Typically initial power spectrum with spectral index n — —2 (which sho uld model CDM like spectra we 
considered in our simulations at small length scales) produces a v alue of .3 which we will be using in our analysis dColombi. Bouchet & Hernguistl 
ll996l : IColombr Bouchet & Schaeffedll996l : IColombi et al.|[l997l) . The VPF cj>{z) and its two-point analog t(z) both exhibit singularity for small 
but negative value of z,, 

<^(z) = <;/>* -a,r(cj,)(z- 2,)""*,; r(2) = r* -b.(2-2,)"'^*"\ (A5) 

For the factorizable model of hierarchical clustering the parameter cj, take s the value —3/2 and a, and b, can be expressed in terms of the nature 
of the generating function GsiT) and its derivatives near the singularity r, teernardeau & Schaeffejl 19921) : 



-| 3/2 

; b, = 



2g's{T,)gn r,) 



-I 1/2 

(A6) 
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As mentioned before the parameter kg which we have introduced in the definition of Qs(t) c an be related to the parameters a and ui appearing in 
the asymptotic expressions of (l){z) tealian & Schaeffeilll989l J Bernardeau & Schaeffeilll992h 

c. = fc„/(fc„ + 2),; a=^^fc^'"="+^ (A7) 

Similarly the parameter Zs which de scribe the behavior of the function 4>{z) near its singula rity can be related to the behavior of Q{t) near Ts 
which is the solution of the equation tealian & Schaeffeilll989l J Bernardeau & SchaeffeJI 19921) 

T* = e^(r.)/e^ (r.), (A8) 

finally we can relate ka to y» by following expression (see Eq. ( IA7I >): 

r. _ 1 {kg + 2f-+^ 

The newly introduced variable x-i, will be useful to define the large-5 tail of the PDF p[5) and the bias 6(5). Different asymptotes in 4>{z) are linked 
with behavior of piS) for various regimes of 5. For very large values of variance i.e. ^2 it is possible to define a scaling function p{5) — h(x) 
which will encode the scaling behavior of PDF, where x plays the role of the scaling variable an d is defined as a; = (1 + 5)7^2. We list different 
ranges of 5 and specify the behavior of p[5) and b[5) in these regimes tealian & Schaeffedl 19891) 

-^/(i-c.) ^ ^ -f, a + 2 ,,,, /a;\i/2 r(a;) /l + (5V'""'^' 



The integral constraints satisfied by scaling function are 5*1 — xh{x)dx = 1 and 5*2 = x'^h{x)dx — 1. These take care of normalization 
of the function p{5). Similarly the normalization constraint over h{5) can be expressed as Cn = xb{x)h{x)dx — 1, which translates into 
J^i dSb{S)p{S) — and f °^ dSSb(S)p(S) — 1. Several numerical studies have been conducted to study the behavior of h(x ) and b{x) for different 
initial conditions see e.g. jColombi. Bouchet & Hemguistl 1 99^ IColombi" Bouchet & SchaeffeJI 1996l ; IColombi et aljl997l) . For very small values 
of 5 the behavior of p{5) is determined by the asymptotic behavior of ^(z) for large values of y, and it is possible to define another scaling function 
g{z) which is completely determined by lj, the scaling parameter can be expressed as z' = {1 + 5)a~^^^^~'^^ However, numerically 

it is much simpler to determine lo from the study of cr{y) compared to the study of g{z). 



i + P(5) = a V 2^^z'a+-)/- ^"Pi'"(T^) J' ^^^^ = "ll"j I— 

(A12) 

To summarize, we can say that the entire behavior of the PDF p{S) is encoded in two different scaling functions, h{x) and g{z') and one can also 
study the scaling properties of b{S) in terms of the scaling variables x and 2 in a very similar way. 



A2.2 The Quasi-linear Regime 

The first departur e from Gaussianity can b e studied analytically using perturbative techniques teernardeaij|l992ll 19941) as well as using steepest 
descent methods l lValageas 2000anli l2002h . These techniques are also inv aluable in modelling the statistics of extreme underdensities. The PDF 
and bias now can be expressed in terms of Gs (t) te ernardeaull 1 9921 . 1 1 994) : 

Gsir) = g{r) -1 = S. (A14) 

The above expression is valid for 5 < Sc where the Sc is the value of 5 which cancels the numerator of the pre-factor of the exponential function 
appearing in the above expression. For S > 5^ the PDF develops an exponential tail whic h is related to the presence of singularity in (j}{z) in a 
very similar way as in the case of its highly non-linear counterpart teemardeaiJI I99II 1994h : 

3a.\/l,-5/2„ „r ,,S \^s\-\^. 1 (1 + 5) 



The quasilinear regime remains relatively better understood than the highly nonlinear regime, as the connectio n to dynamics c an be made using 
analytical schemes. Attempts have been made to extend the perturbative results to the highly nonlinear regime IColombi et alj jl997h . The lower 
order cumulant correlators are related to the moments of the bias function. For the case of top-hat smoothing and specral index rie the lowest order 
of the cumulant conelator is given by teemardeau| [T993): C^f = 68/21 - {n, + 3)/3. 
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APPENDIX B: THE LOGNORMAL DISTRIBUTION 



In addition to the hierarchical ansatz lo/^normal distribution too is a popular analytical model commonly used in cosmology ( iHamiltonlfTgsl ; 
IColes & Jones '1991: Bouche t et ailll993l : | Kofman et al ] | 19941 : iBemardeau & Kofmai]ll995l : IColombilll994h . It appears in the quasi-linear regim e 



jMunshi. Sahni. StarobinskvMl994l) as a natural outcome of gravitational dynamics, under certain simplifying as s umptions i Matarresse et al.ll992h . 



Lognormal distribution is routine ly used to model the statistics of weak lensi ng observables llMunshilbOOOl : IXaniva et al.ll2002h , clustering of 



Lyman-a absorption systems (e.g. iBi & Davidsonl ( Il997l) ) and more recently by jMunshi et alj2011ah for modelling of the tSZ statistics 

To understand the construction of lognormal distribution, we introduce a Gaussian PDF in variable s:p{s) = (27rH)^^''^ exp[— (s — /i)^/2H]. 
With a change of variable s = ln(r) we can write down the PDF of y which is a lognormal distribution p{r) = (27rS)~^^'^ exp[— (ln(r) — 
n)^ /2'E]/r. The extra factor of (1/r) stems from the fact: dr/r = ds. Note that s is positive definite and is often associate with p/{p) = 1 + 5 
which means (r) = 1. The moment generating function for the lognormal interms of the mean fi and the variance S has the following form: 
(r") = exp(Ti/i + n^H/2). This however leads to the fact that if the underlying distribution of r or the density is Gaussian we will have to 
impose the condition: ^ = —3/2.. In our notation above H is the distribution of the underlying Gaussian field. The variance of r defined as 
(r^) — (r)^ = exp(H) — 1 = ^2- So we can write H — ln(l + ^l)- This is the result that was used a bove. The generalisation to two-point or 
bivariate PDF can be done following the same arguments and can be found in jKavo. Taruva. Sutoll200lh . 
In the limit of large separation H12 — >■ we can write down the two point PDF 

PUSi,S2) = Pln{5l)pin{S2)[l + bln{Sl)^Lbln{S2)]; bln{5,) = A, /E . (Bl) 

However, it is simpler to estimate the cumulative or integrated bias associated with objects beyond a certain density threshold Sq. This is defined 
as b{6 > So) = p{5)b{8)d5 / p(S)dS. In the low variance limit ^2 the usual Gaussian result is restored b{S) = The parameters 

A, Ai, S12, H that we have introduced above can be expressed in terms of the two-point (non-linear) correlation function ^^2 = {Si52) and ^2 is 
the volume average of the non-linear two-point correlation function ^12 of the smoothed density field. 

The validity an d limitations of various a spects of the one-point and two-point PDFs have been studied extensively in the literature against 
N-body simulations iBemardeaiJ jl992L IT994}) . It is known that in the weakly nonlinear regime the lognormal distribution is equivalent to the 
hierarchical model with a generating function G{t) = exp(— r). This leads to Sp = p^^^. The loop level corrections can be computed exactly 
for the lognormal distribution, which gives = 3 + and 5*4 = 16 + 15^2 + 6(Cl)^ + It has been shown (see e.g. lKavo. Taruva. Sutol 

that the lognormal distribution very accurately describes the cosmological distribution functions even in the nonlinear regime ^2 < 4 for 
relatively high values of density contrast S < 100. 
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